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The f-electron spectral function of the Falicov-Kimball model is calculated within the dynamical 
mean-field theory using the numerical renormalization group method as the impurity solver. Both 
the Bethe lattice and the hypercubic lattice are considered at half filling. For small U we obtain a 
single-peaked f-electron spectral function, which -for zero temperature- exhibits an algebraic (X-ray) 
singularity (|cj|~ a ) for ui ~ * 0. The characteristic exponent a depends on the Coulomb (Hubbard) 
correlation U. This X-ray singularity cannot be observed when using alternative (Keldysh-based) 
many-body approaches. With increasing U, a decreases and vanishes for sufficiently large U when 
the f-electron spectral function develops a gap and a two-peak structure (metal-insulator transition) . 

PACS numbers: 71.27. +a, 71.20.Eh 



I. INTRODUCTION: 

One of the simplest lattice models for strongly cor- 
related electron systems, the Falicov-Kimball model 
(FKM) 1 , consists of two types of spinless electrons, de- 
localized band (o) electrons, localized /-electrons and 
a local Coulomb (Hubbard) interaction between c- and 
/-electron at the same site. The FKM Hamiltonian reads 

ij i i 

+uJ2 c hflfi ■ (i) 

i 

where Sf denotes the one-particle energy of the /- 
electrons, e c the band center of the c-electrons, the c- 
electron hopping between site i and j and U the Coulomb 
(Hubbard) correlation. This model can be interpreted 
as a simplified version of the one-band Hubbard model 2 
where c and / translates to the label of the different 
spin orientations, and the hopping term for the /-spin 
type is set to zero. The FKM was originally introduced 
as a model for metal- insulator and valence transitions 1 . 
It can also be interpreted as a model for crystallization 
(identifying the "heavy" /-particles with the nuclei, the 
c-particles with the electrons and for U < 0)£. Further- 
more, the FKM is of interest for academic reasons, be- 
cause it is the simplest, non-trivial lattice model for cor- 
related electron systems, for which certain exact results 
are available; for a recent review see Ref. 0. Therefore, 
the FKM may serve as a test case and benchmark for any 
many-body method and approximation, because such a 
method can also be applied to the FKM and comparisons 
of the approximate and the available exact results may 
help to judge the value and possible limitations of such 
many-body techniques. 

Brandt and Mielsch^ solved the FKM exactly in the 
limit of infinite spatial dimensions, d — » oo£, by map- 
ping the problem onto an effective single site subject to 
a complex bath. Thereby, the complex bath field mimics 
the dynamics of the surrounding electrons on the lat- 
tice. Since the local /-electron number is conserved, 



[H, fjfi] = 0, the c-electron problem is solved indepen- 
dently from the dynamics of the /-electrons which are 
static. However, the presence or absence of an /-electron 
on the effective site induces a significant change for 
the conduction electrons: the presence of an /-electron 
switches on an additional scattering potential U for the 
conduction electrons. Brandt and Mielsch^ were able to 
derive an exact, explicit expression for the self-energy of 
the c-electrons S c (z), which is a functional of the local 
c-electron Green function and the (fixed) local /-electron 
occupation numbers. This functional is essentially equiv- 
alent to that obtained in the Hubbard-III-approximation 
of the Hubbard models, i.e. this approximation becomes 
exact for the c-electron self-energy of the FKM in the 
limit d — > oo . 

But despite the conservation of the /-electron occu- 
pation, the /-electrons have a highly non-trivial dynam- 
ics. For U — the f-electron spectral function is simply 
a delta-function, but for finite U the spectral function 
broadens due to the interaction with the fluctuating c- 
electrons. The /-electron Green function Gf(z) of the 
FKM is much more involved and less trivial than the 
c-electron Green function G c , in particular a Hubbard- 
III-type self-energy functional does not become exact 
for the /-electrons. Even for large dimensions, where 
the mapping on an effective impurity problem holds, it 
is not possible to derive an explicit analytical expres- 
sion for the /-electron self-energy functional S/. On 
the contrary, the evaluation of E/ is nearly as compli- 
cated as the self-energy of the full Hubbard model. For 
d — ► oo, the first explicit numerical evaluation of Gf 
(or the /-electron spectral function p f , respectively) was 
performed by Brandt and Urbanekt They started from 
a nonequilibrium, Keldysh based many-body formalism 
using a discretization along the Kadanoff-Baym contour 
and a subtle analytical continuation to obtain the de- 
pendence on real frequency numerically. Only relatively 
few explicit results were presented in Ref. [Tj. Therefore, 
this Keldysh based analysis was repeated very recently 
by Freericks et al£ presenting much more explicit re- 
sults and investigating also the accuracy of the results by 
checking sum rules and the dependence on the numerical 
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discretization and the cutoff in the time integrals. 

On the other hand, it is clear that the problem 
to calculate the /-electron Green function G/(t) = 

— (T(/i(r)//)) must be related to the X-ray threshold 
problem, because the creation of a core hole, with which 
the conduction electrons then interact, is related to the 
creation of an f-electron at a certain time as pointed out 
by Si et al& and Moller et al. lQ As in the X -ray problem 
the non-interacting Fermi sea of c-electrons has to relax 
into an orthogonal ground state through an infinite cas- 
cade of particle-hole excitations after an /-electron has 
been added to the system by fj. However, though the 
relation to the X-ray problem is known, no observations 
of the characteristic X-ray singularities were made in pre- 
vious studies of G/ for the FKM.L& 

In this paper, we apply a different method to the 
FKM, namely the dynamical mean-field theory (DMFT), 
i.e. the (exact) mapping of the d — > oo problem on an 
effective impurity problem^ii, and Wilson's numerical 
renormalization group approach (NRG) 12 as the impu- 
rity solver. Within this DMFT/NRG the typical X-ray 
singularities, in fact, occur for G/ of the FKM. The spec- 
tral function of G/ exhibits a power law \to\~ a divergence 
for the half-filled homogeneous solution at the chemical 
potential when T — > 0. Since the norm of the spectral 
function is finite, the upper limit for a is one: a < 1. 
This new finding has not been previously discussed in 
the existing literature on the /-spectral function of the 
FKM 7, £. The power-law scaling of the /-spectra is re- 
lated to an infra-red problem which -due to the loga- 
rithmic discretization- can properly be resolved within 
the NRG. We will show that the exponent a is propor- 
tional to the imaginary part of the bath field. Since small 
frequencies correspond to a very large time scale, this X- 
ray threshold behavior cannot be seen in the Keldysh- 
contour approach by integrating the equation of motion 
in real time up to a finite cutoff in time2*£. For suffi- 
ciently high temperatures, where the X-ray singularities 
no longer exist, our results are in good agreement with 
Ref . H 

Our investigation can also be motivated as a test or 
benchmark of the different possible (numerical) many- 
body methods, which can be applied to the FKM (and 
thus to correlated electron models, in general). To re- 
solve fine structures on a small frequency scale (around 
the chemical potential) the DMFT/NRG is clearly supe- 
rior to the Keldysh based method. Furthermore, the nu- 
merical efforts and resources needed are smaller by mag- 
nitudes; whereas supercomputers had to be used in Ref|a, 
our results can be obtained on conventional PCs (or lap- 
tops) needing only seconds of computation time. Fi- 
nally, as a byproduct we obtain also the c-electron Green 
function G c and self-energy S c within the DMFT/NRG 
approach, which is known independently from the ex- 
act self-energy functional. Though it is known that the 
NRG is slightly less accurate in the reproduction of high- 
frequency features^ (just because of the logarithmic dis- 
cretization and the resulting lower resolution at higher 



frequencies), very good agreement is obtained also for 
this quantity. 



II. THEORY 

A. Dynamical Mean Field Theory 

The FKM was solved exactly for d — > oo by Brandt and 
Mielsch 5 . Using the locality of the self-energy they 
showed that the local conduction electron Green function 
for the homogeneous phase must have the form 

r loc h] = 1 - n f n f M ) 

c[ > z + ^-£ c - T(z) z + Ji - e c - Y{z) - U> 

since the local /-electron number is conserved: 
[H,fjfi] =0. nf = (fjfi) is the average number of 
/-electrons at site i, and the complex dynamical field 
T(z) contains the dynamics due to the coupling of the ef- 
fective site to the remainder of the system (lattice) . The 
local Green function G c (z) can also be obtained from the 
self-energy E c 



S c (z) = z + f i-T(z)-[G c {z)Y 



via Hilbcrt transformation 



G l c att (z) 



de- 



z + fi — e — S c 



G l c att (z) = G l ° c (z) 



(3) 

(4) 
(5) 



For given /-occupation nf, Eqns. 12151) form a set of 
self-consistency equations from which the dynamical field 
T(z) and the self-energy S c (z) can be calculated. To this 
end, the level position Sf has to be adjusted such that 
the average /-occupation is given by rif . For a fixed to- 
tal particle number n to t per site - n to t — nf + n c - the 
chemical potential is adjusted accordingly. 



B. Dynamical Properties of the Effective Site 

The effective site can also be viewed as an effective 
single impurity Anderson model (SIAM)ii*i5ii£»il 

H eff = (e c - M)c f c+ (e/ - m)/ 1 / 
+ J de ed\d e + Un c hf 

+ [ deVyft^(e)(4c + Jd £a ) . (6) 



with an energy dependent hybridization function A(e) = 
V 2 p eff (e) = 3mr(e — i0 + ) describing the coupling of the 
c-electron to a fictitious bath of "conduction electrons" 
created by d^ with density of states (DOS) p etf (e). The 
/-electron does not interact with the bath field. The 
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hybridization strength V is chosen to be constant and 
defined via 



nV z = 



J deA(e) 



(7) 



and the DOS of the (fictitious) (d) conduction electrons 
is given by p af{ (e) = A(e)/(irV 2 ). 

Using the equation of motion for Fermionic Green func- 
tions, 

z<A\B»(z) = ({A,B})+<g.[A,H]\B^> (z) (8) 

with the commutators [c,H e ff] and [f,H e ff], it is 
straightforward to derive two exact relations for the local 
c- and /-Green function of the effective site, G c (z) =<C 
c|c+> and G f (z) =</|/ f > (z) 

(z-e c -T(z))G c (z) = 1 + U «c/t/| c t» ( z ) (9) 
(z-e f )G f (z) = 1 + U «/ct c |/t» («) (10) 

Parameterizing the Green functions via a self-energy 
S c //, G c = [z — e c - T(z) - E c (z)] _1 and Gf = [z - 
Ef — £/(z)] _1 , yields the exact relations 



£ c (z) = U 



G c (*) 



(11) 
(12) 



with F c (z) =<c/t/| C t> and =</ct C |/t>. The 

problem to be solved is just a special case of a SIAMi& 
with the hybridization matrix element for one spin com- 
ponent set to zero. 



C. The Numerical Renormalization Group (NRG) 

The Hamiltonian © is solved using the NRGi£. The 
core of the NRG approach is a logarithmic energy dis- 
cretization of the effective ci-conduction band around the 
Fermi energy, ui^ — ±D\~ n and A > 1, and a uni- 
tary transformation of the basis d ± onto a basis such 
that the Hamiltonian becomes tridiagonal. The first 
Fermionic operato^ boacr is defined as &o = \ E n ( c L+ + 
d -). Only &o couples directly to the impurity degrees 
of freedom. Eq. © is recasted as a double limit of a 
sequence of dimensionless NRG Hamiltonians: 



n 



lim lim 

A^l+ N^oo 



{d a a-( n -^h n } , 



(13) 



with D A = D(l + A~ 1 )/2, and 



H 



TV 



A" 



+ (V65c + H.c 



N-l N 

Yl A " f ^« { b l+i b n + H.c.} + ^A"*ry„6t6„ 



V is related to the hybridization width V — DaV 

nV 2 = J deA(e) (14) 

while the pre-factor A^ N ~ X '' 2 guarantees that the low- 
lying excitations of TCn are of order one for all N— . 
The energy dependent hybridization function A(w) de- 
termines the coefficients £ n and //,^*-^ by performing 
the unitary transformation from the basis d u ± to b n . 

The NRG is used to calculate the spectral function of 
the Green function for the Fermionic operators A and B 



PA,b{u) = -3to[Ga,b(w - i0 + ) 
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-f3E„ 



-nm - LJ 'mn 



xS(iJ + E n - E rn ) , 



(15) 



where the matrix elements A nm = (n\A\m) are calculated 
using an eigenbasis of the Hamiltonian H, H\n) — E n \n). 
The spectral functions of G c /f and F c /f are obtained 
from their Lehmann representation by broadening the 8- 
function of the discrete spectrum on a logarithmic scale, 
i. e. S(lu -E)-> e- b2 / 4 e-^°^/ E )/ b '> 2 /{^b\E\)^. The 
broadening parameter is usually chosen as 0.5 < b < 1, 
and we used b = 0.8. In order to minimize the broadening 
error, the physical Green functions are calculated indi- 
rectly using the self-energies E c /y(z) given by Eqs. i|TT|) 
and H12|! 18 . The self-energy E c is redundant at this point 
since it was already obtained from the self-consistency 
cycle Eqs. (|2l5p . However, it can be used to quantify 
the discretization error of the NRG, especially at high 
energies^. Since the /-electrons do not hybridize, only a 
spinlcss bath must be considered which makes the NRG 
calculations extremely fast. 



III. RESULTS 
A. Spectral functions for T — > 

We investigate the FKM for the particle-hole sym- 
metric case in the absence of a symmetry break- 
ing using two different densities of states (DOS) for 
the unperturbed conduction electrons on the under- 
lying lattice, namely a semi-circular DOS pa(e) = 
^4 - (e/t*)/(2nt 2 ) (Bethe lattice) and a Gaussian DOS, 
po(e) — exp(— (e/t*)2)/ (i/nt*) (hyper-cubic lattice). 
The energy unit will be t* throughout the paper. At 
half filling, n to t = 1, and for a given set of parameters 
fi = U/2, e c = Ef = 0, we solved the Eqns. 12151) self- 
consistently to determine the bath field T(z). Its imagi- 
nary part A(u>) = 3mr(a; — «0 + ) is used in the algorithm 
described iniS to determine the parameter £ n which en- 
ters the NRG of the effective site of the problem; one has 
rj n = in the particle-hole symmetric case for all n. 
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FIG. 1: Spectral function of the c-Green function of the ef- 
fective site (solid line), the lattice c-Green function, Eqn. fl). 
(dashed line) both calculated using the NRG, in comparison 
with the exact solution obtained from Eqns. 12151 (dotted- 
dashed line) for (a) a Gaussian DOS po(s) and U=l and (b) 
a semi-circular DOS p (e) and U = 2 for T —> 0. NRG pa- 
rameters: A = 1.6, N = 500, b = 0.55. 



As mentioned before, the self-energy S c provided by 
the NRG is a redundant quantity. However, we ob- 
tain G l c att (NRG) from Q as well as G l c oc (NRG) via 
G l ° c (NRG) = [z + n - e c - T(z) - Sc(z)]" 1 (cf. ©) 
which we compare with the exact result from the self- 
consistency equations I|2I5|) . Their spectral functions are 
plotted in figure Q for two different values of U and two 
different DOS pa{e). Fig. Ufa) shows a comparison for 
U = 1 using the Gaussian DOS po while the spectral 
functions for the semi-circular DOS and [7 = 2 are plot- 
ted in nib). The agreement is very good. For \u>\ — > 0, 
the three spectral functions coincide within the numer- 
ical error, a sign of a highly accurate calculation of the 
spectral function for low frequencies due to the equation 
of motion. For large frequencies, the NRG spectral func- 
tions tend to be broadened too much which also yields 
an overestimation of the imaginary part of the self-energy 
S c . Therefore p l f t (NRG) = %mG l c att (NRG)(u;-iO + )/TT 
turns out to be slightly smaller than p l ° c (NRG) for large 
frequencies. This error is controlled by the NRG dis- 
cretization parameter A and can be reduced by reducing 
A while simultaneously increasing the number of states 
N s kept after each NRG iteration. For our calculation, 
we used N s = 500 and A = 1.6. We like to emphasize 
that we only need to add a spin-less chain link in each 
NRG iteration which increases the number of states by a 
factor of 2 contrary to the usual single impurity Ander- 
son or Kondo model where 4 states are needed to describe 
a new chain link which makes the NRG calculations ex- 
tremely fast as mentioned before. 

Two discretization errors determine the deviation of 
the high energy curves: the discretization errors in the 
determination of the NRG hopping parameters £ n by av- 
eraging A (id) over the energy intervals D[A~ n ,A~ n+1 ] 



and the artificial broadening of the 5-function in the 
Lehmann representation of the spectral function l|15l) . It 
is know that the latter is partially compensated using 
the equation of motion^. These systematic errors occur 
in all applications of the NRG as solver of the effective 
single-site problem in the dynamical mean field theory. 
However, these errors do not have any influence on the 
spectral properties for \ui\ — > 0. Here, the NRG furnishes 
an unprecedented resolution due to the logarithmic dis- 
cretization of the energy mesh. This enables us to reveal 
the X-ray properties of the /-spectral function not vis- 
ible in the previously published calculations^, because 
the -in principle exact- method used there requires a dis- 
cretization and time cutoff in the numerical evaluation 
limiting the accessible resolution. 





FIG. 2: /-spectral function for different U for the semi- 
circular Po(b). (a) shows the spectra on a linear- log scale, 
(b) the same data on a log-log scale for positive frequencies. 
For \u)\ — > 0, the spectral function has an algebraic singu- 
larity: Pf(u) cc The inset of (b) displays the lin- 
ear dependence of a on the dimensionless coupling constant 
g — V/o c ff(0) for the metallic case g > 0. With increasing U, 
g decreases monotonically. NRG parameters as in Fig.0 

After we have established the excellent agreement be- 
tween the NRG solution for the FKM and exact result 
for the c-electron spectral function, we will focus on the 
/-spectral function in the remainder of the paper. The /- 
spectral function for the FKM with semi-elliptical DOS 
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is shown in Fig.[3]for T — > 0. Clearly visible in the linear- 
log plot|2ta) is the divergence of the spectral function for 
\lo\ — > 0, while Fig. Efb) reveals its algebraic nature for 
the metallic case (g > 0): Pf(u>) oc |oj|~ q . The /-Green 
function describes the response of the system to the addi- 
tion or removal of an /-electron onto a lattice site. This 
process suddenly changes the Coulomb potential for the 
itinerant conduction electrons. Since the /-electron num- 
ber is conserved in the absence of a i//-hopping matrix 
element, the bath conduction electrons have to relax to 
the new ground state under the presence of the time inde- 
pendent additional scattering potential in a similar way 
as in the X-ray edge problem^. This infinite cascade 
of particle-hole excitations generates a logarithmic slow- 
down which yields an algebraic singularity in the spectral 
function pf{ui). The connection between these two prob- 
lems is mentioned in for instance. But previous 
calculations of the /-spectral functions^ do not reveal 
this feature since the method used there is based on a 
real time integration along Keldysh contours. In order 
to see the algebraic singularities, numerical integrations 
at T = to exponentially large time scales would have 
to be performed, which is not possible in practice also 
because of the possible accumulation of numerical errors. 

The inset of Fig. Elb) shows the dependence of the 
exponent a on the dimensionless coupling constant g — 
Vp eH (0) = 3mr(0 - iO + )/(Vtt) at the chemical poten- 
tial. With increasing U, g is monotonically decreasing. 
At the same time, the exponent a is also reduced as can 
be seen in Fig.^b). We find a linear connection between 
g, which is implicitly U dependent, and the exponent 
a. A relation between U and a was already discussed 
in Ref@. There -for the (unphysical) assumption of a 
Lorentzian unperturbed DOS- a mapping on Nozieres' 
original X-ray problemSLS^ was made. But obviously 
the exponent a obtained within our DMFT/NRG treat- 
ment approaches at a critical interaction strength U c , 
at which the metal-insulator transition sets in, whereas 
no critical U c is obtained in Ref0, probably because only 
a single-impurity (and no selfconsitent DMFT-scheme) is 
considered there. Since the total spectral weight of pf is 
1, the integral over a small interval around zero, on which 
the power law behavior is valid, must be finite, too, i. c. 



dwpf(uj) = 



du>\ 



< 1 



(16) 



Therefore, a < 1 must hold. The exponent a obtained by 
a power-law fit to the spectral function up to |w| ~ 10~ 8 
is consistent with this analytical requirement. 

Figure [3| shows the evolution of the /-spectral func- 
tion Pf(w) for a Gaussian unperturbed DOS po(j-o) as 
function of U for the same NRG parameter as in Fig. [21 
In the chosen units, the square of the effective hybridiza- 
tion strength, V 2 = 0.5 is half the size of the value for V 2 
for the semi-circular DOS. Therefore, the metal-insulator 
transition occurs at a smaller value of U. Setting aside 
this difference and the slightly different shapes, the spec- 
tra are similar. As expected, the spectral functions also 





FIG. 3: /-spectral function for different U for the Gaus- 
sian DOS po(e). (a) shows the spectra on a linear log scale, 
(b) the same data on a log-log scale for positive frequencies. 
For \u)\ — > 0, the spectral function has an algebraic singu- 
larity: p/(tA>) oc The inset of (b) displays the lin- 
ear dependence of a on the dimensionless coupling constant 
g = Vpcft(O) of the effective site. NRG parameters as in Fig.0 



show a power law behavior for \lu\ — * 0. 

We have also checked the basic sum rules pf has 
to fulfill^. Whereas the basic rules / du)pf(u) — 
1,/ du>pf(u>)f{uj) — 0.5 and / dujpf{w)uj — are exactly 
fulfilled, the rula£ / dw p f {uj)uj 2 — C/ 2 /4 is slightly vio- 
lated within our NRG-scheme. We obtain a proportion- 
ality to U 2 but with a larger prefactor, probably as the 
NRG is less accurate at higher frequencies which get a 
higher weight due to the w 2 -factor under the integral. 



B. /-Electron Self-energy 

The imaginary part of — i0 + ) is shown in Fig.0] 
for 4 different values of U and a semi-circular unper- 
turbed DOS. For the metallic case, U — 0.5, 1, 1.8 < U c , 
Qm'Ef(uj~i0 + ) — * for lo — * as expected from the sup- 
pression of scattering processes on conduction electrons. 
However, obviously we observe here again a non-analytic, 
power-law behavior. Therefore, the /-electron spectral 
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FIG. 4: Imaginary part of the /-electron self-energy 5jmS(oj — 
i0 + ) for different values of U and a semi-circular unperturbed 
DOS. The scattering rate increases with increasing U as ex- 
pected from Eqn. 112H . Parameters as in Fig. [5] 



function and self-energy of the FKM provides for an ex- 
ample and paradigm of a system exhibiting non Fermi 
liquid behavior. As it is well known, for a normal Fermi 
liquid the self-energy imaginary part has to vanish ac- 
cording to an w 2 -behavior at the chemical potential. But 
here we have obviously no w 2 -behavior but a non-analytic 
vanishing of QtoE / (ui) . Via G/ = [z — e/ — 
this non-analytic self-energy behavior translates to a di- 
vergence in the spectral function for \u\ — * 0. Above 
a critical value of the Coulomb interaction U, namely 
for U = 2.2 > U c , the solution of the DMFT yields an 
insulating behavior, i.e. a gap in the /-electron spec- 
tral function. Then, the imaginary part of S/(u; — i0 + ) 
has an additional 5(ui) peak, i.e. there — i0 + ) has 

a pole at lu — 0, which leads to a vanishing /-spectral 
function in the energy-gap around the chemical potential. 
This metal insulator transition observed here is similar 
to the Mott-Hubbard transition occuring for the Hub- 
bard model and for the c-electron spectral function of 
the FKM. 




FIG. 5: Evolution of the finite temperature spectral func- 
tion pf(u) for a semi-circular DOS in the metallic phase for 
U = 1.5. We truncated the p(w)-axis in favor of the high 
temperature spectral information. While for high tempera- 
tures, the peak at uj = is rather wide, the development of 
the singularity is clearly visible for the smallest temperature 
plotted here, T = 0.01. NRG parameters as in Fig. Q 




FIG. 6: Evolution of the finite temperature spectral function 
Pf(io) for a semi-circular DOS in the insulating phase for U = 
2.5. While for high temperatures spectral weight is still found 
in the gap between the two sub-bands, it is rapidly vanishing 
for T < 0.1. NRG parameters as in Fig. ID 



C. Finite Temperature 

We used a newly developed finite temperature algo- 
rithm, as described in appendix lAl to calculate the NRG 
spectral functions which enter the Eqns. for the local self- 
energies, Hill and (|12fl . Figure [S] shows the /-spectral 
function for different temperatures in the metallic phase 
for U = 1.5. At high temperatures, the peak at u> = 
is rather broad. Lowering the temperature, the algebraic 
singularity clearly gradually develops. Since we have es- 
tablished its existence in the normal metallic phase in 
the previous section, we truncated the height of the p(uj)- 
axis in favor of the high temperature spectral informa- 
tion. These results agree nicely with those obtained re- 
cently by Freericks et. al£, compare their Fig. 6. Since 
they obtain their spectral function at finite temperature 



by integration in the time domain, their spectral func- 
tion is resolution limited by the largest time ti of the 
numerical integration, uj m i n = 1/t;. For T — ► 0, they 
extrapolate p/(0) to a finite value, which is essentially 
displayed in the inset in Fig. 6 of Ref. |8|. However, our 
work shows that pj(uj = 0,T — * 0) — * oo with a power 
law pf{oS) oc |w|~ Q which could not be seen in Ref. H be- 
cause of the limited frequency resolution inherent to the 
method used. FigureEldisplays the evolution of the finite 
temperature /-spectral function in the insulating phase 
for U = 2.5. While at higher temperature, finite spectral 
weight is found in the vicinity of u> — 0, rather rapidly 
a gap develops for temperatures T < 0.1. The difference 
between the T = 10~ 2 and the T = spectral function 
is marginal for \ui\ > T = 10~ 2 . This behavior is again 
essentially in agreement with results obtained in Ref. 0, 
compare Fig. 7 in that reference. 
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IV. DISCUSSION AND CONCLUSION: 

We calculated the f-electron spectral function pf of 
the FKM at half filling in the d — > oo limit using the 
DMFT/NRG method, pj is strongly temperature de- 
pendent, for small values of U , it has a single peak struc- 
ture with a power law singularity |w|~ Q for lu — > at 
zero temperature T = 0. For larger U > U c , a metal- 
insulator (Mott-Hubbard) transition occurs and pf gets 
a two-peak structure with a gap at T = 0. Therefore, the 
obvious connection of the FKM with the X-ray threshold 
problemSiifl has been explicitcly worked out and demon- 
strated here. The imaginary part of the f-electron self- 
energy vanishes for oj — > for small U < U c , but it does 
not follow the standard Fermi liquid cj 2 -law, but it van- 
ishes non-analytically, i.e. one obtains non Fermi liquid 
behavior—, which leads to the power law singularity in 
Pf. Above the metal- insulator transition, a < and the 
threshhold Eq = f(U — U c ) being a function of U — U c is 
shifted to a finite energy. At finite temperatures, precur- 
sors of the algebraic signularities are found which are cut 
of by thermal fluctuations below \u>\ < T. Furthermore, 
the present investigation demonstrates the strength and 
accuracy of the DMFT/NRG method. It is not only able 
to resolve low frequency or temperature features with un- 
precedented accuracy, it is also very efficient numerically 
requiring only minimal resources and computation time; 
typical run-times have been around 2 minutes on a Pen- 
tium M laptop. 
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APPENDIX A: CALCULATION OF THE FINITE 
TEMPERATURE GREEN FUNCTION USING 
THE NRG 

We used a slight modification of the algorithm for fi- 
nite temperature Green functions by Bulla et. al^i. Our 
algorithm has originally been developed for multi-band 
models. In this case, collecting all matrix elements which 
contribute to pa,b (I15fl at each NRG iteration a nd p atch- 
ing them together as described in Fig. 2 of Ref. |2l| is not 
practical in the two-band model: the number of matrix 
elements exceeds rather rapidly even the memory of mod- 
ern supercomputers for accurate calculations. The patch- 
ing of information from different iterations, however, is a 
very useful concept since (i) it increases the accuracy in 
combination with (ii) producing additional pole positions 
by shifting the eigenenergies slightly through iteration. 
Latter helps to mimic a continuum of states. 

Instead of averaging individual matrix elements for dif- 
ferent iterations, we average different p 1 ^ b (oj) for fixed 
frequency u> but different iterations M. The basic idea is 
to combine the method of Costi et. ali^ with the algo- 
rithm of Bulla et. al. 21 . The eigenstates of iteration M 
represent excitations on an energy scale cum 

lom = D 1 -±^A- M ^ (Al) 

For each iteration M , we define a logarithmic grid inter- 
val Im for 3L + 1 frequencies 



Im — u>m x 



A- 1 ,A- (L - 1)/L ,---,A i = A (J 
I 



-£)/£. 



,A,---,A< 



(A2) 



and evaluate p^ B (u>M^-i) via (|15fl for a fixed iteration 
M. Since the information between even and odd itera- 
tions are collected separately, Im overlaps with the Im—2, 
i.e. I M H Im-i = u M x [l, • • • , A 2 ] and I M H Im-4 = 
lu m x [A, ••■,A 2 ]. Only the p 1 ^ b (lu) evaluated at the 



frequencies I' M = ujm x [1, • • ■ , A] provide accurate infor- 
mation. At each iteration, we add the new data points 
P%B(ui) at Ui(M) = lu m A^- l ^ l , i € [0,3L] to the pre- 
viously obtained set of spectral information by 



- w L 



4 B («iW) 



p AlB (wi(M))- 

^3L - U L 

p% B {wi{M)) for i < L 



U 3 L - Wt 



for i > L 



(A3) 



We do the same for the spectral information at the neg- 



ative frequencies — u>i(M). With this recursion relation, 
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we mimic the patching of the residua since residua far 
away from the frequency ±Wi(M) will only contribute 
marginal through the broadening procedure. It turned 
out that L = 3 — 7 is a good choice, and for this work, we 
used L = 6. Since the spectral information B (uji(N)) 
for i = 0, . . . , L — 1 is not very reliable for the last NRG 
iteration N, we drop these L frequencies. 



The major advantage of this new algorithm is that any 
temperature dependent spectral function can be calcu- 
lated on the fly There is no need for storage of informa- 
tion of previous iterations other than the spectral func- 
tion pa,b(u) itself. We also use two different broadening 
functions 



-&V4 e -(l°g("/E)A>) 2 /(^,|£|) f or \ E \ > LwT 



depending whether the pole position E is larger than a 
multiple of the temperature T, E > L W T or E < L W T . 
L w is a constant of the order of 1 as suggested in the 
literature^. The parameter Lt controls the width of 
the Lorentzian on an energy scale of T. Here, we used 



L w = 0.3 and Lt = 0.1. Since the iteration M also 
represents a characteristic temperature scaled, we only 
include spectral information up to iterations^ for which 
lom > Lu er T , where we chose Li ter = 0.1. 
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